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The dynamics of two coupled piece-wise linear one-dimensional monostable maps is investigated. 
The single map is associated with Poincare section of the FitzHugh-Nagumo neuron model. It is 
found that a diffusive coupling leads to the appearance of chaotic attractor. The attractor exists in 
an invariant region of phase space bounded by the manifolds of the saddle fixed point and the saddle 
periodic point. The oscillations from the chaotic attractor have a spike-burst shape with anti-phase 
synchronized spiking. 
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Spiking-bursting activity is a common feature 
of the temporal organization in many neural firing 
patterns [1]. Bursting activity means that clus- 
ters of spikes (of action potentials) occur more or 
less rhythmically and separated by phases of qui- 
escence. Spiking-bursting dynamics can be reg- 
ular or chaotic depending on the concentration 
of neuromodulators, currents and other control 
parameters. For example, many of the talamo- 
cortical neurons from central patterns generators 
generate chaotic spiking-bursting dynamics. In 
order to understand the emergence of chaotic os- 
cillations in a neurones network, we use a two 
variables FitzHugh-Nagumo model of the mem- 
brane potential of an isolated neural cell. For 
appropriate values of the parameters, the system 
possesses a stable fixed point a focus surrounded 
by a separatrix such that for large enough per- 
turbations, the system responds by an excita- 
tion pulse which also decays to the stable focus. 
Hence, introducing a section by a half-line, we 
obtain a one-dimensional piece-wise smooth map 
with one discontinuity point representing the ex- 
citation threshold. Then we consider a difference 
coupling between two such maps. In the linear 
approximation of the maps and for suitable values 
of the parameters, this two-dimensional system is 
locally hyperbolic having discontinuity lines and 
a strange attractor . The interesting fact is that 
these lines define regions in the phase space corre- 
sponding to the thresholds of excitability of each 
one of these neurones, or none of them. The sym- 
bolic dynamics represents then the bursting ac- 
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tivity of the neurones. It is then possible to ap- 
ply the approximation tools of statistical analysis 
of time series generated by the system, namely 
through Perron-Frobenius Operator, in order to 
study and interpret some important features of 
neural networks, like synchronization and time 
correlations. 



I. INTRODUCTION 

Many neural systems composed of a number of inter- 
acting neurons exhibit self-sustained oscillatory behavior 
leading to the formation of various space-time patterns. 
Such patterns are believed to play a key role in signal 
and information processing functions of the brain 
One of the fundamental problems is the understanding of 
possible dynamic mechanisms of such patterns to appea r 
and to evolve in time and space [H IE 01 IS EU Il3.ll3| . 
There are two basic phenomena of emergence of oscil- 
lations (regular or chaotic) in neuron assemblies. The 
first one is obvious and deals with the presence of lo- 
cal infra-cellular oscillations. Being coupled such units 
are capable of various oscillatory patterns, clustering and 
synchronization. The second one, found in recent theo- 
retical and experimental studies, concerns the possibil- 
ity of oscillations in assemblies of non-oscillatory cells 
IS ■ The oscillations may also appear in coupled 
non-identical cells for sufficiently strong coupling. The 
assembly is characterized by an oscillatory "average" cell 
dynamics which makes the non-oscillatory cells oscillat- 
ing. Another studies have reported that coupling even 
identical excitable cells can modify the dynamics to form 
oscillatory attractors co-existing with a stable fixed point 
• Such attractors are characterized by anti-phase spik- 
ing. The effects of electrical coupling between neural cells 
also include the appearance of bursting in two coupled 
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pacemaker cells, modification of burst period in coupled 
bursters, synchronization and chaos [8lMITol [T^ . [l3Lll^| . 

A model approach, in order to display the dynami- 
cal origin of neural oscillations, is to use a simplified 
behavior-based description of the system. For this pur- 
poses nonlinear maps could be helpful as they can pro- 
vide an appropriate qualitative description of complex 
dynamic processes including chaotic behavior in lower di- 
mensional systems 0, EH Ell ■ In this paper we study the 
system of two coupled piece-wise linear one-dimensional 
maps. The single map is derived from the dynamics 
of an isolated neural cell modelled by the FitzHugh- 
Nagumo excitable system 0, 0, The FitzHugh- 

Nagumo neuron model can be derived from Hodgkin- 
Huxlcy conductance-based model, for some parameter 
values, when we take into account the difference of ki- 
netics between the potential dependent gating variables 
0- The first variable describes the evolution of neuron 
membrane potential, the second mimics the dynamics of 
outward ionic currents. Then, the model can describe the 
salient features of neuron dynamics including the action 
potential generation, excitability and excitation thresh- 
old. The map has one globally stable fixed point and a 
discontinuity corresponding to the excitation threshold of 
the cell. We shall show how linear diffusive coupling be- 
tween the two maps leads to the appearance of chaotic os- 
cillations with anti-phase spiking. A number of studies of 
coupled chaotic maps have shown that anti-phase chaotic 
oscillations may appear when the synchronization mani- 
fold looses transverse stability and the system evolves to 
the off-diagonal attractors [11 El El El HI . 

The paper is organized as follows. In Sect. II we show 
how the dynamics of excitable FitzHugh-Nagumo model 
can be described by a piece-wise continuous point map 
and introduce two-dimensional map modeling a pair of 
coupled cells. In Sect. Ill we analyze the dynamics of the 
map. We numerically show the emergence of strange at- 
tractor in an invariant domain defined by invariant mani- 
folds of saddle fixed point and saddle periodic orbit. Sect. 
IV describes the statistical characteristics of the attractor 
set and the emergent chaotic oscillations with anti-phase 
spiking. Section V contains a brief discussion of the re- 
sults. 



II. POINT MAP DESCRIPTION OF THE 
EXCITABLE DYNAMICS 

To replicate the excitable dynamics of an isolated neu- 
ral cell one can use the FitzHugh-Nagumo-like model. It 
can be taken in the following form 

u = u — m 3 /3 — v, 

v = e(k(u) — v — I), ^ 

The u-variable describes the evolution of the membrane 
potential of the neuron, v describes the dynamics of the 
outward ionic currents (the recovery variable) 3]. The 




Fig. 1 

FIG. 1: Qualitative picture of the phase plane of Eqs. S 
illustrates the transversal half-line. 

function k is taken piece-wise linear, k(u) = k\u, if 
u < 0, and fc(tt) = k2U, if u > 0. The parameters 
k\ and control the shape and the location of the v- 
nullcline, hence the dynamics of the recovery term. The 
parameter e defines the time scale of excitation pulse 
and the parameter / is a constant stimulus. The ex- 
citable behavior of Eqs. Q is illustrated in Fig. 1. 
Appropriate values of the parameters provide the exis- 
tence of three fixed points 0\ [u (V) , v^), 2 (^ 2 \w^) 
and Oz(u( 3 \ v^). The points 0\ and O3 are stable and 
unstable foci, respectively, the point O2 is a saddle with 
the incoming separatrix defining the excitation thresh- 
old. Then, if a perturbation of the rest state, 0\, is large 
enough, i.e. lies beyond the separatrix, the system re- 
sponds with an excitation pulse, otherwise it decays to 
the stable rest point 0\ (Fig. 1). 

In order to describe the dynamics of the unit by 
a point map, we introduce the transversal half-line L : 
{v = v^\u > u^} (Fig. 1). It is found that the return 
map of the flow given by Eqs. at the section L defines 
a map, $ : L — > L, for all points excluding the one, Lq, 
corresponding to the intersection of the incoming saddle 
separatrix with the half- line L. This point never returns 
to the cross-section. Consequently, we also must exclude 
all the pre-images of the point Lq, = (&~ 1 ) m Lo, m = 
1,2, . . . , 00 . Then, the Poincare map $ can be written 
as 

x(n+l) = ${x(n)), (2) 

n = 1,2,..., 00, 

with x(n) accounting for (u — it^ )-coordinates of the 
points at the Poincare section. The map is invertible 
and defined in the interval [0, 00) excluding the points 
Lfe . The shape of the curve calculated numerically 
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FIG. 2: The curve &(x) obtained from Eqs. Q at the section 
half-line L. Parameter values: e = 0.5, ki = 0.5, k-z = 2,1 = 
0.221. 



III. DYNAMICS OF THE MAP 



A. General properties of map 



Since F(x) is given by the piece-wise linear function, 
then the differential Df of the map / is a constant ma- 
trix, with eigenvalues 



Mi 



2d, /z 2 



a 



(5) 



Hence, Lyapunov exponents are defined by A; = 
ln|/Zj|,i = 1,2, for a arbitrary trajectory of the map, 
except zero Lebesgue measure set of initial conditions. 

Let us consider the system (0 in the parameter region 



D = {a, d, a : < a < a/b, 



±4^ <d< 1, 



>0} 



For convenience, we shall fixed a = 1. This region 
D is chosen in order to have Lyapunov exponents Ai > 
0, A 2 < 0, hence any trajectory, {f l xo, xo £ B}, where 
B is a set of discontinuous points of the map /, is heavily 
dependent on initial conditions, and to have the modulus 
of Jacobian of the map / less than one, in D. All these 
features allow us to speak of chaotic behavior. 

We shall study the trajectories of the map / missing 

oo 

the discontinuity set B = [J f^ 1 ^, where 



is shown in Fig. 2. It is given by a piece- wise continuous 
curve with one stable fixed point. Then, the dynamics 
of the map is trivial. All its trajectories represent se- 
quences of points monotonically approaching to the fixed 
point. The discontinuity point plays the role of the ex- 
citation threshold, i.e. the neuron exhibits spiking if the 
map evolves above this point. Note, that such map ob- 
tained from gives an "average" description of the cell 
dynamics. 

Let us approximate the piece-wise continuous curve 
by a piece- wise linear function and define the map 
for all points (— oo, oo) excluding the points Lfc. Then, for 
two cells we introduce difference coupling term (electri- 
cal coupling) and consider the following two-dimensional 
map /: R 2 — ► R 2 



x 1 (n + 1) = F(xi(n)) + d(x 2 (n) — x\{n)), 
x 2 (n + 1) = F(x 2 {n)) + d(xi(n) - x 2 (n)), 



(3) 



where the variables x\^ 2 refer to the dynamics of the two 
cells, d is the coupling coefficient, the function F(x) is 
taken in the form 



F(x) 



ax, if x < a; 

ax + a{b — a) if x > a, 



(4) 
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= {xi,x 2 : {xi = a} [J{x 2 = a}} 



(6) 



For convenience, let us change variables, according to 
the eigenvectors of Df 

yi (rc) = x 2 {n) - xi{n), y 2 (n) = x 2 (n) + xi(n) (7) 

Then, the discontinuity lines become 

Ti = {2/1,2/2 : 2/2 -2/1 = 2a} 
T2 = {2/1,2/2 : 2/2 + 2/1 = 2a} 

The lines Ti and T 2 divide the phase plane (2/1,2/2) of 
the map / into four regions 



G\ = {2/1,2/2 
G 2 = {2/1,2/2 
G3 = {2/1,2/2 
Ga = {2/1,2/2 



2/2 + 2/1 < 2a, 2/2 - 2/1 < 2a} 

2/2 + 2/1 > 2a, 2/2 - 2/1 < 2a} 

2/2 + 2/1 < 2a, 2/2 - 2/1 > 2a} 

2/2 + 2/1 > 2a, 2/2 - 2/1 > 2a} 



In each region Gi the map / is continuous and has the 
form 



with 



< a < 1. 



2/1 (n+ 1) = Mi2/i(«-) + b\ 
y 2 (n+l) = fi 2 y 2 {n) + b 2 , 



0, if (2/1,2/2) E Gi,G 4 
a(b-a), if (2/1,2/2) € G 2 
-a(b - a), if (2/1,2/2) 6 G 3 



(8) 
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0, if (2/1,2/2) G Gi 
o 2 = <( a(& - a), if (2/1, 2/2) G G 2 , G 3 
2a(6-a), z/ (2/1,2/2) G G 4 

The map / has one hyperbolic fixed point, O(0, 0). In 
region G\ its invariant curves coincide with the coordi- 
nate axes. Analyzing map / in the regions G2, G4 we 
find that it has a hyperbolic periodic orbit, Q, of period 
two with coordinates 

a(b-a) _ a{b-a)_ 

2/1 = —j 1 7=2/l(0), 2/2 = —j = 2MG), 



2d - a - 1 



1 - a 



if the parameters satisfy the inequality 



(9) 



j j . , t j 2a-a 2 (a + b) 

d<d p with d p = — \ ■ 10 

2[2a — a(a + ojj 

Stable and unstable invariant manifolds of the orbit Q 



are 



n --,ni - J „ ■ Vi=Vi{Q)' l f (2/1,2/2) G G 2 
W(0)-<! 2/1,2/2. tfl = _ yi (g),i / ( !ft , !tt ) e G, 



^"(Q) = { 2/1,2/2: 2/ = 2/ 2 (<9), i/ (2/1,2/2) GG 2 ,G 2 -} 
It follows from @ that 

lim 2/i(<3) = +°°- 

Then periodic orbit Q appears from infinity in the phase 
plane (z/i, y 2 ) for the parameters belonging to the curve 



1 - , , , 1 + a a 

N_ 1 = <a.d: d= , < a < - 

1 2 b 



in the parameter plane (a,d). Note, that on this curve 
one of the multiplier hits the bifurcation value, [i\ = 
— 1. Then fixed point O changes its stability becoming 
of saddle type in D. 



Let us introduce region II in the phase plane defined 
by 

n = {1/1, 2/2 : -2/1 (0) < 2/1 < 2/i (0), < y < y 2 {Q)} 

Let us find the conditions for region II to be an in- 
variant region of map /. Since the system (jSJ is sym- 
metric with respect to j/ 2 -axis it is sufficient to con- 
sider only 2/1 > region. Assuming (2/1(0), 2/2(0)) £ G±, 
where G+ = {II \ B} f| d f|{2/i > 0}, let us obtain 
the conditions on the parameter values for which 
(2/i(l), 2/2(1)) € II. In region Gf we find that 



yi (l) = -(2d-a) yi (0) 
2/2(1) = 0:2/2(0) 



(13) 



It follows from JED that 2/2(1)) G II if the follow- 

ing inequalities 

Vl[Q) < 2/1(0) 2/2(0) <^ (14) 



2d - a 



2d -a 



are satisfied. Using the symmetry of / with respect to 
2/1, we can restrict the above condition to the right part 
of II, and since 



(15) 



max y 1 (0)<y 1 (M), 
max y 2 (0) < y 2 (Q) 

(s/i(0),3/ 2 (0))eG+ 

the inequalities (|14f> are satisfied if 



f^->2/l(M), ^> 2/2(0) (16) 
2a — a a 

The first inequalities in l|l(jfl impose the following condi- 
tion on d: 



d < dh with dh 



l + 2a+A/l 



2a(b-a) 



(17) 



B. Invariant region and chaotic attractor of the 
map / 

Let us consider the location of invariant curves of the 
points O and periodic orbit Q in the phase plane with 
respect to the lines Ti, T 2 . The unstable invariant curve 
(the separatix) Wi(Q) intersects the line T 2 at point K 
with coordinates 

2a — a(b + a) , „ 

2/1 = — ^ -, 2/2 = 2/20, 11 

1 — a 

and invariant curve W"(0) at the point M with coordi- 
nates 



2/1 = 2a, y 2 = 0, 



(12) 



Figure 3 illustrates the location of the invariant curves of 
fixed point O and orbit Q in the phase plane (2/1,2/2)- 



Thus, under the condition (|1T|1 /(G^) C II. Similarly 
we find that the image of G^ = {II \ B} f| G 2 f| {t/i > 0} 
by / is also included in II, for the parameter in D. 

Finally, II is invariant region of map / in the parameter 
region 

D mv = D f]{d < d p } f]{d < d h } 

Figure 4 shows region D inv in the parameter plane 
(a, d). In this plane the boundary of Di nv consists of 
three components: 

N_i = {a,d : d= ±±«, < a < f } 
N b = {a,d: a= f, ^ < d < d h ) 
N h ={a,d:d = d h ,0<a<^}. 

The line N_i corresponds to the appearance of periodic 
orbit Q and changing stability of the fixed point O. The 
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Fig. 3 

FIG. 3: Mutual location of invariant manifolds of the fixed 
point O and periodic orbit Q, and discontinuity lines ]?i,2 on 
the phase plane of map /. 
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FIG. 5: Chaotic attractor A on the phase plane (1/1,2/2) 
rameter values: b = 4.95, a = 0.2, d = 0.74. 
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FIG. 4: Parameter region Di nv on the parameter plane (a, d). 



line TVf, is the boundary of the monostability of the uncou- 
pled, d = 0, maps ©■ The points of the curve Nh corre- 
sponds to the "tangency" of the separatricies IF"(0) and 
Wf(Q), and for these parameter values f(M) £ W{(Q~). 

Therefore, if the parameter values of system © be- 
longs to region Di nv , then invariant (absorbing) region 
II and /(II \JS) C II exists in the phase plane. Conse- 
quently, II contains strange (chaotic) attractor A of the 
map /. Figure 5 illustrates possible structure of the at- 
tractor A in the phase plane. 



To characterize the complexity of the chaotic attractor 
A we calculated numerically its fractal dimension d f (A) . 
It appears that, indeed, df(A) takes non-integer values 
greater than 1 (Fig. 6). Note that the box dimension of 
the set increases with the increase of coupling coefficient 
d. Corresponding estimate of the dimension df using 
Lyapunov exponents is shown by dashed curve. 



C. Chaotic oscillations and attractor. 

Figure 7 (a) illustrates time evolution of the variables 
x\ (n) , X2 (n) corresponding to chaotic attractor A. Let 
us characterize the oscillations in terms of original model 
describing neuron excitability. Note that the neuron 
excitation threshold accounted by saddle separatrix (Fig. 
1) corresponds to the discontinuity lines T\, T2 in the map 
description. Thus, if a map trajectory jumps above these 
lines then we can refer this event as a neuron spike, if it 
evolves below Ti, T2, the neuron is not excited . In such 
a way the map oscillations can be described with binary 
variables z;\ 



Zi(n) 



0, if Xi(n) < a; 

1, if Xi(n) > a. 



The evolution of z^-variablcs is shown in Fig. 7 (b). It 
appears that map oscillation represent sequence of bursts 
containing anti-phase spiking. Note that characteristic 
time scales of oscillations (burst duration and inter-burst 
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FIG. 6: Fractal dimension df(A) of the chaotic attractor A 
calculated by box counting method versus coupling coefficient 
d. Parameter values: b — 4.95, a = 0.2. 



(a) Fig- 7 




s5 1 iio 



FIG. 7: (a) Chaotic oscillations of the map ©. 'xi,2(n) are 
shown by solid and dashed curves, respectively, (b) Spike- 
burst behavior of binary variables z\^{n). Parameter values: 
b = 4.95, a = 0.2, d = 0.74. 



period) depend coupling d. For smaller d one can find 
longer lasting subthreshold period and shorter burst pe- 
riod. Then, the map describes a kind of chaotic spike- 
burst behavior typical for many neural systems. 



IV. STATIONARY PROBABILITY 
DISTRIBUTION ON THE CHAOTIC 
ATTRACTOR 

The global behavior of the chaotic dynamics is de- 
scribed by the probability of finding the trajectory in 
any given region of the attractor. It can be visualized by 
a distribution of a cloud of points, each moving under the 
deterministic mapping. A stationary distribution may be 
reached by the system after long term evolution. To ob- 
tain this distribution, we start from a histogram of clouds 
of points in each small "box" of the phase space. These 
points are mapped by / and a new cloud is thus obtained. 
Thus to the initial probability density po correspond a 
new density p\. The operator which maps po — > p\ is 
called the Perron- Frobenius operator: p\ = Ppo defined 
by: 



Ppo(yi,y2)dyidy 2 



Po(yi,y2)dyidy 2 (18) 



for any region A of the phase space. The stationary den- 
sity p is an eigenfunction of P corresponding to the eigen- 
value 1. The existence and uniqueness of the stationary 
distribution has been the object of many works [1,112 ■ It 
is nevertheless difficult to obtain analytical exact expres- 
sion of p(yi, 2/2), apart from the restricted case of Markov 
maps. We shall use an approximating algorithm inspired 
from the method of Ding and Zhou [24J . 

The phase space of the coupled system is divided into 
n x m identical rectangles A^. We consider an initial 
probability density p Q (2/1,2/2) that is constant on each 
A,.: 



Po(Y) = J2 



Oi{0) 
m(A,-) 



XAi (Y) 



(19) 



where Y — (2/1,2/2), rn(Ai) is the Lebesgue measure of 
Aj and cii(0) is the probability of finding a phase point 
in A,- 



h, XAi is the characteristic function of Aj: 



if 
if 



Y e A t 
Yd A, 



It is clear that p\ has no reason to be of the "Coarse- 
grained" form l|19fl . Smoothing this density by integrat- 
ing it on each Aj we obtain by using the definition |^ 



aj (l) = J A .Pp (Y)dY 

= f f - H A ]} Po(Y)dY 



= E^f^XA^dY 

m (f 1 (Aj) n 

m(Ai) 



(20) 



^(A ^ ) 



Thus, the transition probability from Aj to Aj is given 
by the stochastic matrix: 



Pid = rn(f 1 (A j )/A i 



m(Aj 



)) (21) 
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The matrix p = (pij) is an approximate value of the 
Perron-Frobenius operator P. It can be proved that p 
converges to P as Aj — > 24]. Eq. lj2"n)l can be 
written under the matrix form: 

a(t+l) = a(t)p, t = 0,... (22) 

where a(t) is the row vector (oi(t), . . . ,apf(t)) with N = 
n x m. Thus, an approximate stationary probability is 
the row eigenvector v = (vx(n), . . . , uj\r(n)) such that: 

u = u • p, (23) 

We first calculate the matrix elements of p and then 
we compute the eigenvector v. To calculate pij we use 
1)21 [I. In each rectangle Aj we put fc x £ points uniformly 
distributes - then we compute the number of point (Y £ 
Aj) such that f(Y) G Aj, and we divide this number by 
the number of points in Aj. This provides ptj. 

We used this method to approximate the stationary 
distribution for a — 0, 2 and d = 0, 75 (Fig. 8) showing 
a density distributed over some region, which is to be 
compared with the Fig. 5 where we have used only one 
trajectory. It is to be noted that the last figure shows 
the attractor generated by one trajectory whereas Fig. 8 
shows the approximate density of the attractor. 

V. SYNCHRONIZATION. 

The above method allows to obtain the statistical dis- 
tribution of synchronized spikes. Recall that the vari- 
able x\ is spiking if and only if Y € and 22 is 
spiking if and only if Y S G 2 , and no one is spiking if 
Y G Gf U Gi which we denote G\. Recall also that the 
region G4 corresponds to x\ and X2 are simultaneously 
spiking, which lies outside the invariant region. More- 
over, because Gx, G 2 are disjoint, x% and X2 cannot be 
simultaneously spiking . Thus we shall consider a parti- 
tion of the invariant domain II corresponding to one of 
the three possible states of spiking S(Y): 

S(Y) = x u if Y eGz 
S(Y) = x 2s if YeG+ 
S(Y) = if Yed, 

where xi s means that the only spiking variable is Xi\ and 
means that no one is spiking. This partition induces 
symbolic dynamics, that we shall study. We calculate 
the stochastic matrix 7Tjj from any one of these states 
to any other, where i,j are the states x ls , a^s, 0, for two 
values of the coupling parameter d = 0.65, 0.7 (see Fig. 
9). This table shows that we have only the following 
probable transitions in one step: 

X 2s ► Xi s > X 2s 

and < — > 0. 

Less probable transitions are: 

xu — ► 0, x 2s — ► 



Xls 




X 2s 

but these transition probabilities depend monotonically 
on the coupling d. Thus the transitions x\ s — > X2 S and 
X2 S — ► Xu (i.e. successive spikes) decreases slightly with 
coupling. So, strengthening the coupling in this range 
decreases the bursting probability. 

We can understand these occurrences in considering 
the intersections of the image of each region by / with 
these same regions (see Fig. 10). The most probable 
transition xi s — > Xk s , {k, 1} S {1, 2} 2 and — ► corre- 
spond respectively to f(Gf) f| Gf and f(d) f| G x . The 
less probable transitions xi s — > and — ► xi s cor- 
respond respectively to f(G2±)f)Gi and /(Gi)P)G^. 
The positive Lebesgue measure of these intersected re- 
gions give the transition probabilities of each of these 
occurrences. Some of the transition probabilities given 
by the calculation of the Lebesgue measures of the inter- 
sected regions are shown in Fig. 11, with their depen- 
dance over the coupling constant d, which can be com- 
pared to those obtained numerically in the table Fig. 9. 
All the calculations are made in appendix A, where we 
also show that ,when the parameters are in Dj nt ,,a neuron 
can not fire twice, i.e. G 2 f] f(G 2 ) = 0. 

We can now check the memory effects in studying the 
conditional probability of a state given the past. If we 
denote by So, S-i, S—2 ■ • ■ the values of the variable S at 
time, 0, —1, —2 etc, we like to know if conditional prob- 
ability Prob(S \S- 2 ,S-i) is equal to Prob(S Q \ S-x) or 
not. In the first case there is no memory of the process. 
In order to compute the above conditional probability 
we have to find Prob(So = io | S-2 = i-2,S-i = *-i) 
which implies (3) 3 = 27 calculations. But, neglecting 
rare events like S-2 = 0, S-i = X2 S having 0, 0422 proba- 
bility for d — 0, 75, we have only 4 significant transitions. 

Prob(S = x 2s I 5_2 = x 2s , S-i = xi s ) = Pi 
Prob(S = xi s I S-2 = xu, S-i = x 2s ) = Pi 
Prob(S = I S-2 = x 2s , S-i = xis) = P3 
Prob(S = I S-2 = xu, S-i = x 2s ) = Pi 

The first one is the probability of having X2 S if it is 
preceded by a sequence of x 2s and a;i s and for the others 
the definition is the same. Thus we obtain the following 
table (Fig.12). The fact that the probability Pi = 0,879 
is not equal to Prob(So — X2 S | S-i = xi s ) = 0,9303 
means that the system has acquired a memory, the chain 
is not a Markov chain. The memory of a chain is given 
by the smallest positive integer m a such that: 

Prob(S Q I S- m , S- m +i, ■ ■ -,S-i) = Prob(S | S , _ m0r .., s _ 1 ) 

for any m > m , and any value of (So, ■ ■ ■ S- m ). 

It is also remarkable that the probability of having 
if it is preceded by the sequence (x 2s , xi s ) is 0, 1209 the 
double of the probability of having preceded by xi s 
simply 0, 0696. This explain the frequent appearance of 



FIG. 8: Approximative stationary distribution of the map /.Parameter values: b — 4.95, a = 0.2, d = 0.74. 
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FIG. 9: The 3*3 stochastic matrices of transition probabilities 
for three values of the coupling d. 
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FIG. 10: Image by / of the partitions of the invariant region 

n. 



sequence of spikes. How could we estimate the memory 
of the process depends on the full calculation of all con- 
ditional probabilities for all past. 

Another very interesting aspect is the dependence of 
this memory on the coupling coefficient as shown in Fig. 
13. It is clear that the memory effects increase with cou- 
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FIG. 12: The conditional probabilities of spiking and non- 
spiking after two successive spikes respectively. 
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FIG. 11: Probability transition from a spike to another, and 
from rest to spike, versus coupling constant d\ b — 4.95, a — 
0.2. 



FIG. 13: The effect of coupling on the memory effect. Pa- 
rameter values: b = 4.95, a = 0.2. 



pling (the distance between the two curves increases with 
coupling). 



VI. CONCLUSION 

We have investigated chaotic dynamics of two cou- 
pled maps each formally derived from FitzHugh-Nagumo 
model of neuron excitability. Being uncoupled such maps 
are trivial with only one stable fixed point corresponding 
to neuron rest state. The excitation threshold is given by 
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the discontinuity point. Wc have introduced the formal 
definition of spike if the map evolves above this threshold. 
The linear "diffusive" coupling term introduced for the 
maps can be treated as a kind of electrical inter-neuron 
coupling that mimics an "integral" coupling current be- 
tween the two cells. 

We have shown that increasing the coupling coeffi- 
cient above certain threshold leads to the appearance of 
chaotic attractor. The attractor appears in the invari- 
ant region of the phase space confined by the invariant 
manifolds of saddle periodic orbit. This region attracts 
trajectories from outside and does not contain any stable 
trajectory inside. The oscillations emerging with such 
an attractor have a spike-burst shape with synchronous 
bursts with anti-phase spiking. Using probabilistic de- 
scription we have found the probabilities to find a spike 
in different conditions. It has been also found that the 
system acquires a memory (in contrast with Markov's 
processes). 



Acknoledgments 



we shall show that (/(Gi) f| G*) \J(f(Gi) fl G 2 ) ^ 9, 
which means that a neuron at rest can fire. That is to 
say 

y 1 (/ 1 (M-))>y 1 (M+), 

-2a(a - 2d) > 2a [AZ > 

However, the last inequality is always true when wc 
are in D inv . 

Next, we compute the image of G^~: 

h{M+) = {a(-a + b) + 2a(a - 2d),a(-a + &)}, 
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f (D+\ — r ajb-a) a(b-a) 
J2W ) — l a j_i_2d' 1-q Si 

f2(S + ) = {^- d ,a(b-a)}, 

(A3) 

We first shall see that f{G 2 ) is completely included 
in Gi |J G2 and excluded of G 2 , which implies that one 
neuron can not spike twice. For the simplicity of proof, 
we shall use indirect calculations. We are going to show 
that: 
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APPENDIX A: IMAGE OF THE INVARIANT 
REGION AFTER ONE ITERATION. 

In order to understand the oscillatory behavior it is 
interesting to look at the image of the invariant region 
after one iteration (see Fig. 10). First of all, we have 
to determine the images of the regions of the invariant 
region, defined by the discontinuities, namely Gi, Gj 
and G^"(e.q. to G^). 

Gi is a polygon having for vertices P + , P~ , M + and 
M"; Gt, Q~, P~, S- and M"; and G^ , P+, Q+, 
S + and M~ (when M+ < S + or equivalent, which is 
always satisfied for d < dh)- Let us compute the images 
of all these vertex, according to the definition of / in 
each region. We have to keep in mind that some of these 
points belong to the lines of discontinuity; in this case, 
the computation is taken in the sense of a limit starting 
from the interior of each region. 

First, we compute the image of Gi obtained by the 
image of its vertices: 



/i(M+) = {2o(a-2d),0}, 
/i(M-) = {-2a(a-2d),0}, 

f 1 (P+) = r {a(b+a)-2a)(a~2d) a 2 (a-b) j 

f _ ( (a(j»+a)-2a)(a-2d) a 2 (a-b) 1 

Jl\ r ) — t a -l ' a-1 J' 



(Al) 



Vi(h(M+)) < 0, 
yMK)) > 0, 

MM*)) > vMK))) 

y2(MK))-y 1 (f 2 (K)<2a 



(A4) 



The first three inequalities allow us to say that f 2 {M + ) is 
under and on the left of f 2 {K)\ but, as P + is on the seg- 
ment [K M+], f 2 (P + ) is on the segment [f 2 (K] f 2 (M+)}, 
and so it is under and on the left of /2 (if). The last in- 
equality show that f 2 (K) (and by consequence /2(P + ))is 
not in G2 ■ 

First, we verify the first inequality 



tfi(/ 2 (M+)) < 0, 

a(b - a) + 2a{a - 2d) < 0, 

a(b + a) < 4ad, 



(A5) 



However in D; 



a(b + a) < a + aa, 
2a + 2aa < 4ad, 



(A6) 



but a + aa < 2a + 2aa is always satisfied, as a > 0, 
a > 0. By consequence, 

yi(/ 2 (M+)) < 0. 

The second and the third inequalities are straightfor- 
ward as (b — a) > 0, a > and a > 0. 
Finally, we check the last inequality: 

V2{h{K)) - yi(f 2 (K) = a{b - a) + 2aa - a(b - a), 
- Vi(f2(K) = -2ab <0<2a, 

(A7) 

and so a neuron can not fire twice. 
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Finally, we compute G 2 : 

h{M~) = {-a(-a + b) - 2a(a - 2d), a(-a + b)}, 

r (p-\ _ f afr(2a-l-2rf)-a(a+2rf(g-2)) a(b-a) -| 
J3\ > ~ I 1— a ' 1-a J' 

f fn-l — r MM <*(&-a) i 

J3W ; — 1 a4-l-2d' 1-a J' 

/ 3 (S-) = {-^E^,«(6-a)}, 

(A8) 



By symmetry, the above results for G 2 hold for G 2 ■ 



The computation of the Lebesgue measure of the re- 
gions of intersection involve only computation of polygo- 
nal area which it is straightforward with the coordinates 
of the vertices. 
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